##############################################################################
# File-Name: 07_bayes_factor.r
# Purpose: Produces bayes factor analysis
# Machine: macOS Ventura 13.5.1
# R version 4.3.1 (2023-06-16)
##############################################################################

### Required packages --------------------------
pacman::p_load(tidyverse, scales,  knitr, srvyr, extrafont, 
               rebus, broom, tidyr, lubridate, here, ggdist, ggtext, ggalt, janitor, 
               broom, tidyr,  stargazer, janitor, estimatr, list, marginaleffects, posterior, 
               logspline, effectsize, bayestestR, rstanarm)

# Data wrangling ----------------------------------------------------------

# open
d <- readRDS(here("data","all_processed_data.rds"))


# variables for all primary utcomes----------------------------------------------------------------
misinfo_outcomes = c("false_false_sum",  "true_true_sum", 
                     "false_news_exp", "true_news_e")

list_politems_items <- list("pol_all_zcore", 
                            "distance_candidates",
                            "affective_pol_diff_all", 
                            "social_pol_agg", 
                            "mean_abs_policy_polarization")

list_wb_items <- list("swb_zcore", # zscore index
                      "q_well_being_happy_numeric"  ,
                      "q_well_being_depressed_numeric",
                      "q_well_being_anxious_numeric",
                      "q_well_being_isolated_numeric", 
                      "q_well_being_satisfied_numeric" )

formula = "~ exp"


# Model ------------------------------------------------------------------
# Setting this seed here so that replication is the same. 
# Unfortunately I did not set a seed in the conditionally accepted version. 
# The findings are substantively, even though there is some numerical variation in the simulation
set.seed(1310)


bayes_factor_model <- function(dv, data){
    
    # build model.expr  
    model.expr = as.formula(paste(dv, " ~ exp"))
    
    # fit a bayesian model
    bayesfit <- stan_glm(
        formula = model.expr, 
        data = data,
        prior = normal(0, .5, autoscale = T), 
        prior_intercept = normal(.5, .5, autoscale = FALSE)
    )
    
    # Compute the Bayes factors using the `bayesfactor_parameters` function
    BF <- bayesfactor_parameters(bayesfit, null = 0, direction = "two-sided")
    
    # Compute Bayes factors in the direction of the alternative for all model parameters
    BF$Bayesfactor <- (exp(BF$log_BF))
    
    return(BF)
    
}

# estimate for all mode
bayes_res <- map(c(misinfo_outcomes,list_politems_items, 
                   list_wb_items), ~ 
                     bayes_factor_model(.x, data=d)  %>%
                     as_tibble() %>%
                     mutate(dv=.x))


# Combine everything
bayes_tidy <- bayes_res %>%
    bind_rows() %>%
    filter(Parameter=="exptreatment") %>%
    slice(1:4, 5, 10) %>%
    mutate(labels=c("False News Accuracy", "True News Accuracy", "False News Exposure", 
                    "True News Exposure", "Polarization Index", "Subjective Well-Being Index")) %>%
    mutate(labels=fct_relevel(labels,c( "False News Exposure", 
                                        "True News Exposure", 
                                        "False News Accuracy", "True News Accuracy",
                                        "Polarization Index", "Subjective Well-Being Index")),
           Qualit=NA) %>% arrange(labels)

library(kableExtra)
kable(
    bayes_tidy %>%
        select(labels,Bayesfactor, Qualit), 
    caption = "Bayes Factors in the direction of the alternative hypothesis (BF10)",
    format = "latex", booktabs = T, escape = F, linesep = "",
    row.names = F, label = "bayes") %>%
    kable_styling(full_width = F, latex_options = c("HOLD_position", "scale_down")) %>%
    column_spec(3, width = "7cm") %>%
    footnote(general = "",
             threeparttable = T) %>%
  save_kable(file ="output/sm_tab20.tex")
